#################################################################################################
#################################################################################################

###          R code to replicate findings presented in Supplementary materials                ###
###     (Appendix B, sections B3-B4, and Appendix D: Tables B6-B9 and Figures B4 and D1)      ###
###        to "The Index of Emancipative Values: Measurement Model Misspefication"            ###
                          

#################################################################################################
#################################################################################################



#################################################################################################

###                          I. LOADING DATA AND REQUIRED R PACKAGES                          ###

#################################################################################################


library(foreign)
library(MplusAutomation)


#WVS<- read.spss("C:\\Users\\lssi\\Desktop\\pro_choice values #invariance\\WVS_Longitudinal_1981_2014_spss_v2015_04_18.sav",
#                to.data.frame=T, 
#                 use.value.labels=T)

## WVS <- subset(WVS, select=c(S002, S003, S017, S019,
## F118, F120, F121, E012))

## Names for relevant variables
## S002: The wave
## S003: The country
## S017: The original weight
## S019: The original weight equilibrated to an homogeneous N of 1500 for each 
## Choice: F118, F120, F121  (Homosexuality, Abortion, Divorce)
## E012: Willingness to fight for one's own country
 
WVS<- read.spss("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\WVS_choice.sav",
                to.data.frame=T, 
                 use.value.labels=T)



#################################################################################################

###                      II. RECODING AND RENAMING RELEVANT VARIABLES                         ###

#################################################################################################


WVS$F118 <- as.numeric(WVS$F118)
WVS$F120 <- as.numeric(WVS$F120)
WVS$F121 <- as.numeric(WVS$F121)

### Create a new variable "Cultural zone" (as defined by Welzel [2013]) 

WVS$Zone[WVS$S003 == "Algeria" | 
WVS$S003 == "Iran" | 
WVS$S003 == "Iraq" |
WVS$S003 == "Libya" | 
WVS$S003 == "Saudi Arabia" |
WVS$S003 ==  "Egypt" |
WVS$S003 ==  "Jordan" | 
WVS$S003 == "Morocco" |
WVS$S003 == "Turkey"|
WVS$S003 == "Palestine" |
WVS$S003 == "Libya"|
WVS$S003 ==  "Yemen" |
WVS$S003 ==  "Lebanon" |
WVS$S003 == "Tunisia" |
WVS$S003 == "Bahrain" |  
WVS$S003 == "Qatar"|            
WVS$S003 == "Kuwait"] <- "Islamic East"# , 16 countries

WVS$Zone[WVS$S003 == "Bangladesh" | 
WVS$S003 =="Indonesia"|
WVS$S003 =="Pakistan" | 
WVS$S003 =="India"|
WVS$S003 =="Malaysia" |
WVS$S003 == "Philippines"|
WVS$S003 =="Singapore"|
WVS$S003 =="Thailand"] <- "Indic East"#, 8

WVS$Zone[WVS$S003 == "Viet Nam"|
WVS$S003 =="China"|
WVS$S003 =="Hong Kong"|
WVS$S003 =="South Korea"|
WVS$S003 =="Taiwan"|
WVS$S003 =="Japan"] <- "Sinic East"  #, 6

WVS$Zone[WVS$S003 == "Azerbaijan"|
WVS$S003 =="Albania"|
WVS$S003 =="Armenia"|
WVS$S003 =="Belarus"|
WVS$S003 =="Bosnia"|
WVS$S003 =="Georgia"|
WVS$S003 =="Kyrgyzstan"|
WVS$S003 =="Russia"|
WVS$S003 =="Serbia"|
WVS$S003 =="Kazakhstan"|
WVS$S003 =="Uzbekistan"|
WVS$S003 =="Ukraine"|
WVS$S003 =="Macedonia"|
WVS$S003 =="Moldova"|
WVS$S003 =="Romania"|
WVS$S003 =="Bulgaria"|
WVS$S003 =="Serbia and Montenegro"|
WVS$S003 =="Montenegro"] <- "Orthodox East" #, 18

WVS$Zone[WVS$S003 == "Cyprus"|
WVS$S003 =="Greece"|
WVS$S003 =="Israel"|
WVS$S003 =="Andorra"|
WVS$S003 =="Austria"|
WVS$S003 =="Belgium"|
WVS$S003 =="France"|
WVS$S003 =="Ireland"|
WVS$S003 =="Italy"|
WVS$S003 =="Luxemburg"|
WVS$S003 =="Malta"|
WVS$S003 =="Portugal"|
WVS$S003 =="Spain"] <- "Old West" #, 13

WVS$Zone[WVS$S003 == "Denmark"|
WVS$S003 =="Finland"|
WVS$S003 =="Germany"|
WVS$S003 =="Iceland"|
WVS$S003 =="Netherlands"|
WVS$S003 =="Norway"|
WVS$S003 =="Sweden"|
WVS$S003 =="Switzerland"|
WVS$S003 =="Great Britain"] <- "Reformed West" #, 9

WVS$Zone[WVS$S003 == "Australia"|
WVS$S003 =="Canada"|
WVS$S003 =="New Zealand"|
WVS$S003 =="United States"] <- "New West" #,4

WVS$Zone[WVS$S003 == "Croatia" |
WVS$S003 =="Latvia"|
WVS$S003 =="Lithuania"|
WVS$S003 =="Estonia"|
WVS$S003 =="Czech Rep."|
WVS$S003 =="Hungary"|
WVS$S003 =="Poland"|
WVS$S003 =="Slovakia"|
WVS$S003 =="Slovenia"] <- "Returned West"  #, 9

WVS$Zone[WVS$S003 == "Guatemala"|
WVS$S003 =="Venezuela"|
WVS$S003 =="Ecuador"|
WVS$S003 =="Colombia"|
WVS$S003 =="Mexico"|
WVS$S003 =="Peru"|
WVS$S003 =="Brazil"|
WVS$S003 =="Chile"|
WVS$S003 == "Puerto Rico"|
WVS$S003 =="Dominican Rep."|
WVS$S003 =="El Salvador"|
WVS$S003 =="Argentina"|
WVS$S003 =="Trinidad and Tobago"|
WVS$S003 =="Uruguay"] <- "Latin America" #, 14 countris

WVS$Zone[WVS$S003 == "Burkina Faso"|
WVS$S003 =="Ghana"|
WVS$S003 =="Ethiopia"|
WVS$S003 =="Nigeria"|
WVS$S003 =="Rwanda"|
WVS$S003 =="Tanzania"|
WVS$S003 =="Uganda"|
WVS$S003 =="Zimbabwe"|
WVS$S003 =="Mali"|
WVS$S003 =="Zambia"|
WVS$S003 =="South Africa"] <- "sub-Saharan Africa"  #10 #,11 coun tries and total 108




#################################################################################################

###                   III. PREPARING DATA FOR MPLUS INVARIANCE ANALYSIS                       ###


#################################################################################################


### WAVE 1 data preparation ###

# identify countries where one or more pro-choice items were not asked
WVS1 <- WVS[WVS$S002 == "1981-1984",]
WVS1$S003 <- droplevels(WVS1$S003)
table(WVS1$S003, WVS1$F118) # None
table(WVS1$S003, WVS1$F120) # None
table(WVS1$S003, WVS1$F121) # None

# create MPLUS data file
prepareMplusData(WVS[WVS$S002 == "1981-1984",], 
"WVS_1.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

# Returns country codes (to be inserted into MPLUS input file)
data1 <- read.table("WVS_1.dat")
table(data1$V2)
# 9   10   56   71   82   87  107  148


### WAVE 2 data preparation ###

WVS2 <- WVS[WVS$S002 == "1989-1993",]
WVS2$S003 <- droplevels(WVS2$S003)
table(WVS2$S003, WVS2$F118) # None
table(WVS2$S003, WVS2$F120) # None
table(WVS2$S003, WVS2$F121) # None


prepareMplusData(WVS[WVS$S002 == "1989-1993",], 
"WVS_2.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

data2 <- read.table("WVS_2.dat")
table(data2$V2)
# 9   22   28   35   36   46   73  
# 82   87  107  121  129  136  144 
# 148  150  155  163 


### WAVE 3 data preparation ###

WVS3 <- WVS[WVS$S002 == "1994-1998",]
WVS3$S003 <- droplevels(WVS3$S003)
table(WVS3$S003, WVS3$F118) # Bangladesh, Pakistan, Turkey
table(WVS3$S003, WVS3$F120) # Bangladesh, Pakistan, Turkey
table(WVS3$S003, WVS3$F121) # Pakistan, Turkey


prepareMplusData(WVS[WVS$S002 == "1994-1998"& WVS$S003 != "Bangladesh"  & 
WVS$S003 != "Pakistan" & WVS$S003 != "Turkey",], 
"WVS_3.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

data3 <- read.table("WVS_3.dat")
table(data3$V2)
# 2    8    9   10   14   24   28   35   36  
# 37   38   42   46   49   51   55   56   58   61
# 71   73   82   87   93   97  107  110  111  118  
# 121  122  127  128  129  133  135  136  140 
# 144  146  148  150  154  155  166  167  169 
# 171  174  176  193 


### WAVE 4 data preparation ###

WVS4 <- WVS[WVS$S002 == "1999-2004",]
WVS4$S003 <- droplevels(WVS4$S003)
table(WVS4$S003, WVS4$F118) # Iraq, Morocco, Turkey
table(WVS4$S003, WVS4$F120) # Turkey
table(WVS4$S003, WVS4$F121) # Turkey


prepareMplusData(WVS[WVS$S002 == "1999-2004" & WVS$S003 != "Iraq"  & 
WVS$S003 != "Morocco" & WVS$S003 != "Turkey" ,], 
"WVS_4.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

data4 <- read.table("WVS_4.dat")
table(data4$V2)
#  2    3    9   13   20   31   35   36  
# 73   74   75   78   82   84   87   89 
# 107  110  111  121  123  127  128  133 
# 138  140  143  145  148  149  150  165 
# 167  168  170  171  176 

### WAVE 4 Reduced 

prepareMplusData(WVS[WVS$S002 == "1999-2004" & WVS$S003 != "Iraq"  & 
WVS$S003 != "Morocco" & WVS$S003 != "Turkey" &
WVS$S003 != "Algeria" & WVS$S003 != "Bangladesh"  & 
WVS$S003 != "Pakistan" & WVS$S003 != "Saudi Arabia" ,], 
"WVS_4_red.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

data4red <- read.table("WVS_4_red.dat")
table(data4red$V2)
#  2    9   20   31   35   36   73   74 
#  75   78   82   84   87   89  107  110 
# 111  121  127  128  133  140  143  145 
# 148  149  150  165  167  168  170  171  176 


### WAVE 5 data preparation ###

WVS5 <- WVS[WVS$S002 == "2005-2009",]
WVS5$S003 <- droplevels(WVS5$S003)
table(WVS5$S003, WVS5$F118) # Iraq, Morocco, Peru, Egypt
table(WVS5$S003, WVS5$F120) # Peru, Egypt
table(WVS5$S003, WVS5$F121) # Peru


prepareMplusData(WVS[WVS$S002 == "2005-2009"& WVS$S003 != "Iraq"  & 
WVS$S003 != "Morocco" & WVS$S003 != "Peru" & WVS$S003 != "Egypt",], 
"WVS_5.dat", keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

data5 <- read.table("WVS_5.dat")
table(data5$V2)
#5    9   10   22   24   31   35   36   37
#38   44   53   56   57   58   61   62   65
#70   71   73   74   75   79   82   84   87 
#101  102  107  110  117 118  122  129  135  
#136  137  145  146  148  150  154  155  158  
#160 163  166  169  171  173  174  178  179 


### WAVE 6 data preparation ###

WVS6 <- WVS[WVS$S002 == "2010-2014",]
WVS6$S003 <- droplevels(WVS6$S003)
table(WVS6$S003, WVS6$F118) # Kuwait and Egypt
table(WVS6$S003, WVS6$F120) # Egypt
table(WVS6$S003, WVS6$F121) # None

prepareMplusData(WVS[WVS$S002 == "2010-2014"& 
WVS$S003 != "Kuwait"  & WVS$S003 != "Egypt",], 
"WVS_6.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

data6 <- read.table("WVS_6.dat")
table(data6$V2)
#   3    8    9   10   12   14   22   28   
#  35   36   37   38   44   50   55   58 
#  60   61   62   70   73   76   82   83
#  84   87   89   91   95  101  107  112 
# 117  118  121  123  127  128  129  134
# 135  136  137  143  146  148  149  150 
# 154  158  160  162  163  166  171  174  175  177 


### WAVE 6 reduced

prepareMplusData(WVS[WVS$S002 == "2010-2014"& 
WVS$S003 != "Kuwait"  & WVS$S003 != "Egypt" & WVS$S003 != "Bahrain" &
WVS$S003 != "Palestine"& WVS$S003 != "Jordan"& WVS$S003 != "Lebanon" & 
WVS$S003 != "Morocco"& WVS$S003 != "Pakistan",],  
"WVS_6red.dat",
keepCols = c("S002", "S003", "S019", "Zone", 
"F118","F120", "F121"))

data6red <- read.table("WVS_6red.dat")
table(data6red$V2)
#3    8    9   10   14   22   28   35  
#36   37   38   44   50   55   58   61 
#62   70   73   76   82   83   87   89  
#95  101  107  117  118  121  127  128 
#129  134  135  136  137  143  146  148 
#149  150  154  158  160  162  163  166 
#171  174  175  177 



### Run MPLUS models (MPLUS imput files were created separately)
runModels("C:/Users/lssi/Desktop/pro_choice values invariance/Appendix 2 MPLUS/Pro-Choice Invariance",
recursive = T, showOutput = T)




#################################################################################################
 
###                                      IV. FIGURES                                          ###

#################################################################################################

opar <- par()

###     FIGURE B4: Relationship between the raw country mean scores for "choice" and the      ###
###            Bayesian country mean scores in WVS Waves 3 to 6.                              ###


# Compute Welzel's version of the "choice" index
WVS$Homosexuality <-(as.numeric(WVS$F118)-1)/9
WVS$Abortion <-(as.numeric(WVS$F120)-1)/9
WVS$Divorce <-(as.numeric(WVS$F121)-1)/9
WVS$prochoice <- (WVS$Homosexuality + WVS$Abortion + WVS$Divorce)/3


# compute raw country means on pro-choice values for WVS waves 3 to 6.

means.agg3 <- aggregate(WVS[WVS$S002 == "1994-1998"& WVS$S003 != "Bangladesh"  & 
WVS$S003 != "Pakistan" & WVS$S003 != "Turkey",]$prochoice,
 by = list(WVS[WVS$S002 == "1994-1998"& WVS$S003 != "Bangladesh"  & 
WVS$S003 != "Pakistan" & WVS$S003 != "Turkey",]$S003), FUN = mean, na.rm = T)

means.agg4 <- aggregate(WVS[WVS$S002 == "1999-2004" & WVS$S003 != "Iraq"  & 
WVS$S003 != "Morocco" & WVS$S003 != "Turkey",]$prochoice,
 by = list(WVS[WVS$S002 == "1999-2004" & WVS$S003 != "Iraq"  & 
WVS$S003 != "Morocco" & WVS$S003 != "Turkey",]$S003), FUN = mean, na.rm = T)

means.agg5 <- aggregate(WVS[WVS$S002 == "2005-2009"& WVS$S003 != "Iraq"  & 
WVS$S003 != "Morocco" & WVS$S003 != "Peru" & WVS$S003 != "Egypt",]$prochoice,
 by = list(WVS[WVS$S002 == "2005-2009"& WVS$S003 != "Iraq"  & 
WVS$S003 != "Morocco" & WVS$S003 != "Peru" & WVS$S003 != "Egypt",]$S003), FUN = mean, na.rm = T)

means.agg6 <- aggregate(WVS[WVS$S002 == "2010-2014"& 
WVS$S003 != "Kuwait"  & WVS$S003 != "Egypt",]$prochoice,
 by = list(WVS[WVS$S002 == "2010-2014"& 
WVS$S003 != "Kuwait"  & WVS$S003 != "Egypt",]$S003), FUN = mean, na.rm = T)


# Extract Bayesian means from Mplus output files
modelout3 <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 3.out")
modelout3p <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 3_partial.out")
modelout4 <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 4.out")
modelout4p <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 4 partial.out")
modelout5 <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 5.out")
modelout5p <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 5 partial.out")
modelout6 <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 6.out")
modelout6p <- readModels("C:\\Users\\lssi\\Desktop\\pro_choice values invariance\\Appendix 2 MPLUS\\Pro-Choice Invariance\\wave 6 partial.out")


means.app3 <- as.numeric(modelout3$parameters$stdyx.standardized[modelout3$parameters$stdyx.standardized$paramHeader== "Means",]$est)
means.app3p <- as.numeric(modelout3p$parameters$stdyx.standardized[modelout3p$parameters$stdyx.standardized$paramHeader== "Means",]$est)
means.app4 <- as.numeric(modelout4$parameters$stdyx.standardized[modelout4$parameters$stdyx.standardized$paramHeader== "Means",]$est)
means.app4p <- as.numeric(modelout4p$parameters$stdyx.standardized[modelout4p$parameters$stdyx.standardized$paramHeader== "Means",]$est)
means.app5 <- as.numeric(modelout5$parameters$stdyx.standardized[modelout5$parameters$stdyx.standardized$paramHeader== "Means",]$est)
means.app5p <- as.numeric(modelout5p$parameters$stdyx.standardized[modelout5p$parameters$stdyx.standardized$paramHeader== "Means",]$est)
means.app6 <- as.numeric(modelout6$parameters$stdyx.standardized[modelout6$parameters$stdyx.standardized$paramHeader== "Means",]$est)
means.app6p <- as.numeric(modelout6p$parameters$stdyx.standardized[modelout6p$parameters$stdyx.standardized$paramHeader== "Means",]$est)


### Create Figure B4

tiff('Figure B4.tiff', units = "cm", width = 24, height = 14, res = 600)

par(mfrow= c(2,4))

cor(means.app3, means.agg3$x) # 0.9177296
plot(means.app3, means.agg3$x, pch = 19, cex = 0.6,
main = 'Wave 3', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg3$x ~ means.app3))
text(-.57, 0.67, expression(paste("Pearson's ", rho, "= 0.918")),
cex = 1)

cor(means.app3p, means.agg3$x) # 0.9164132
plot(means.app3p, means.agg3$x, pch = 19, cex = 0.6,
main = 'Wave 3 (partial)', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg3$x ~ means.app3p))
text(-.58, 0.67, expression(paste("Pearson's ", rho, "= 0.916")),
cex = 1)

cor(means.app4, means.agg4$x) # 0.816759
plot(means.app4, means.agg4$x, pch = 19, cex = 0.6,
main = 'Wave 4', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg4$x ~ means.app4))
text(-1.85, 0.5, expression(paste("Pearson's ", rho, "= 0.817")),
cex = 1)

cor(means.app4p, means.agg4$x) # 0.8176903
plot(means.app4p, means.agg4$x, pch = 19, cex = 0.6,
main = 'Wave 4 (partial)', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg4$x ~ means.app4p))
text(-1.72, 0.5, expression(paste("Pearson's ", rho, "= 0.818")),
cex = 1)


cor(means.app5, means.agg5$x) # 0.9354984
plot(means.app5, means.agg5$x, pch = 19, cex = 0.6,
main = 'Wave 5', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg5$x ~ means.app5))
text(-1.08, 0.75, expression(paste("Pearson's ", rho, "= 0.935")),
cex = 1)

cor(means.app5p, means.agg5$x) #0.9342704
plot(means.app5p, means.agg5$x, pch = 19, cex = 0.6,
main = 'Wave 5 (partial)', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg5$x ~ means.app5p))
text(-1.02, 0.75, expression(paste("Pearson's ", rho, "= 0.934")),
cex = 1)

cor(means.app6, means.agg6$x)  # 0.9154726
plot(means.app6, means.agg6$x, pch = 19, cex = 0.6,
main = 'Wave 6', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg6$x ~ means.app6))
text(-1, 0.75, expression(paste("Pearson's ", rho, "= 0.915")),
cex = 1)

cor(means.app6p, means.agg6$x) # 0.9046571
plot(means.app6p, means.agg6$x, pch = 19, cex = 0.6,
main = 'Wave 6 (partial)', xlab = 'Latent Means Bayes',
ylab = 'Sum Scores')
abline(lm(means.agg6$x ~ means.app6))
text(-1.2, 0.75, expression(paste("Pearson's ", rho, "= 0.905")),
cex = 1)

dev.off()



#################################################################################

###              Figure D1: Nation-level association between                  ###
###               pro-choice values and willingness to fight                  ###

will.agg6 <- aggregate(as.numeric(WVS[WVS$S002 == "2010-2014"& 
WVS$S003 != "Kuwait"  & WVS$S003 != "Egypt",]$E012) - 1,
 by = list(WVS[WVS$S002 == "2010-2014"& 
WVS$S003 != "Kuwait"  & WVS$S003 != "Egypt",]$S003), FUN = mean, na.rm = T)

cor(means.app6, will.agg6$x)   #  -0.4431004
cor(means.app6p, will.agg6$x) #  -0.4433526
cor(means.agg6$x, will.agg6$x) #  -0.4993515

tiff('Figure D1.tiff', units = "cm", width = 24, height = 13, res = 600)

par(mfrow=c(1,2))

plot(means.app6, will.agg6$x, pch = 19, cex = 0.9,
xlab = 'Latent Means on Pro-Choice Values (Bayes)',
ylab = 'Mean Score on Willingness to Fight')
abline(lm(will.agg6$x ~ means.app6))
text(0.9, 0.95, expression(paste("Pearson's ", rho, "= -0.443")),
cex = 1)
text(0.25, 0.3, "Japan", cex = 0.9)
text(1.45, 0.79, "Sweden", cex = 0.9)


plot(means.agg6$x, will.agg6$x, pch = 19, cex = 0.9,
xlab = 'Raw Means on Pro-Choice Values',
ylab = 'Mean Score on Willingness to Fight')
abline(lm(will.agg6$x ~ means.agg6$x))
text(0.65, 0.95, expression(paste("Pearson's ", rho, "= -0.499")),
cex = 1)
text(0.42, 0.3, "Japan", cex = 0.9)
text(0.75, 0.79, "Sweden", cex = 0.9)

dev.off()

#q()
